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Abstract. We study a previously largely unexplored branch of homogeneous and isotropic 
background solutions in ghost-free massive bigravity with consistent double matter coupling. 
For a certain family of parameters we find ‘self-inflated’ FLRW cosmologies, i.e. solutions with an 
accelerated early-time period during the radiation-dominated era. In addition, these solutions 
also display an accelerated late-time period closely mimicking GR with a cosmological constant. 
Interestingly, within this family, the particular case of ( 3 \ = @3 = 0 gives bouncing cosmologies, 
where there is an infinite contracting past, a non-zero minimum value of the scale factor at the 
bounce, and an infinite expanding future. 


Keywords: Cosmology, massive gravity, bigravity, modified gravity, bouncing universes 


Contents 


1 Introduction 1 

2 Bigravity model with double matter coupling 3 

3 Cosmological background 5 

3.1 Branch I 7 

3.2 Branch II 8 

4 Branch I 10 

4.1 Positive case: Nb > 0 10 

4.2 Null case: Nb = 0 14 

4.3 Negative case: Nb < 0 16 

4.4 Radiation and dust 19 

5 Summary, discussion and conclusions 20 

A Branch I: Friedmann equation for effective metric 22 

B Matter-dominated universe 25 


1 Introduction 

Our present cosmological standard model - ACDM and with it General Relativity (GR) - 
still faces several problems. These range from the technical un-naturalness of the cosmological 
constant A to requiring the presence of further unknown ‘dark’ components, namely dark energy 
(if not fully accounted for by A) and dark matter, making up nearly 95% of the total effective 
matter content of the universe [1], if GR is indeed the correct theory of gravity on all scales. 
Massive bigravity, as proposed in [2, 3], is a natural extension to massive gravity [4-6], and a 
promising alternative to GR, which may help to answer some of the problems mentioned above. 
This model is an interacting bimetric theory containing more physical degrees of freedom (DoF) 
than GR, as it propagates two spin-2 fields, corresponding to one massive (5 DoF) and one 
massless (2 DoF) graviton, in contrast to the single massless graviton of GR. 

Massive bigravity passes several of the immediate consistency checks any theory of gravity 
needs to pass. It agrees with solar system constraints on the presence of a fifth force [7— 
9] (the massive graviton would propagate such an extra force) via the Vainshtein mechanism 
[10]. On cosmological scales, it also gives rise to viable homogeneous and isotropic (FLRW) 
background solutions [11-14]. The specific massive bigravity interactions (between both metrics) 
were carefully chosen to avoid an extra propagating scalar DoF, namely a Boulware-Deser ghost 
[2, 4-6, 15-18], causing the theory to be unstable. However, the absence of the Boulware-Deser 
ghost does not guarantee the model to be instability free, as some of the DoF of the gravitons 
may still behave as ghosts, tachyons or create gradient instabilities. In fact, this is generically 
the case, as has been shown at the level of perturbations in previous cosmological studies (see 
below). 

The cosmology of ghost-free massive bigravity has been extensively studied when matter 
is coupled to one of the metrics only (single matter coupling) [11-13, 19-42], In this setting, 
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two main branches of solutions were identified. Even though at the level of the background 
both branches could lead to viable cosmologies, at the level of linear perturbations the so-called 
expanding branch (also known as finite branch) was found to have gradient instabilities on the 
scalar sector, leading to an exponential growth of these perturbations [ 22 ], and in turn breaking 
the validity of perturbation theory at early times 1 . On the other hand, the so-called bouncing 
branch (also known as infinite branch) was found to have the helicity -0 mode of the massive 
graviton behaving as a ghost, and has tachyon instabilities in the tensor sector [33, 34, 38-41] 2 . 

More general consistent theories of metrics/spin-2 fields beyond massive (bi-)gravity with 
a single matter coupling have also been explored recently. New kinetic interactions were investi¬ 
gated in [45-51], generalisations of the potential interactions of massive bigravity to N multiple 
metrics in [52-60], and new couplings to matter in [14, 31, 56, 61-81]. These matter couplings 
allow matter to couple to both metrics and we will therefore refer to them as ‘double matter 
couplings’. Generic such couplings re-introduce the Boulware-Deser ghost at an unacceptably 
low scale, however the specific couplings of [70, 71] stand out in that they are consistent ghost- 
free double matter couplings. Note that this statement is taken in an effective field theory sense. 
They are the only known non-derivative couplings that do not introduce a ghost up to the A 3 
strong coupling scale (for a discussion of derivative couplings see [82]), so massive bigravity with 
such a double matter coupling can be safely considered as an effective field theory with a cutoff 
at A 3 or above. A ghost is ‘present’ at larger energies, which are beyond the regime of validity 
of the theory (in fact the ghost scale may set the theory’s cutoff), but the ghost is of course not 
at all present in the low-energy effective theory below the cutoff. In fact it has recently been 
argued that the couplings of [70, 71] are the unique such couplings for a cutoff > A 3 [83-85]. 

In the context of this double coupling, some homogeneous and isotropic cosmological so¬ 
lutions were studied in [14], where viable (background) evolutions were found. However, at the 
level of linear perturbations, tachyonic, gradient, and ghost instabilities were found for these 
solutions [ 86 , 87] in the tensor, vector and scalar sector, respectively. In the light of these 
results, in this paper we analyse previously unexplored homogeneous and isotropic background 
solutions with the consistent double coupling mentioned above. Specifically, we find new self- 
inflated cosmologies, by which we mean solutions with an accelerated early-time period during 
the radiation-dominated era. Depending on the choice of parameters, this inflating period could 
happen within the regime of validity of massive bigravity, and thus give an interesting alterna¬ 
tive to standard cosmic inflation, where, contrary to GR, no unknown inflationary field would 
need to be introduced. However, even in the case where the inflationary early universe phase is 
outside of the regime of validity of the theory, the new solutions we find also give new consistent 
background cosmologies for massive bigravity at late times, closely mimicking GR with a cos¬ 
mological constant. In addition, all the self-inflated solutions found have a minimum non-zero 
value of the scale factor (and associated finite energy densities) and thus avoid any physical 
Big Bang singularity. The particular case f3\ = ^3 = 0 gives a bouncing universe, i.e. where 
there is an infinite contracting past, followed by a bounce with an infinite expanding future. A 

1 Some papers have suggested ways to overcome this issue [43, 44], though. 

2 General relativity with a cosmological constant is of course a particular limit of ghost-free bigravity and 
massive gravity theories and as such one would naively expect that a GR-like and hence instability-free evolution 
can be recovered by suppressing any non-GR operators, i.e. by establishing a hierarchy between coupling constants 
in the theory. Recent studies have confirmed that this is indeed the case. For instance, in a study of the expanding 
branch in [43], it was found that for appropriate choices of parameters, the scalar exponential instability could 
be moved outside the regime of validity of the theory (in an effective field theory sense), or equivalently to ‘very 
early times’. When curing the instabilities of the theory in this way, the price to pay is of course suppressing any 
non-GR signatures. 
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detailed analysis on the stability of linear perturbations around these backgrounds is beyond the 
scope of this paper, but we argue that some improvements are expected compared to previous 
solutions found in [86, 87]. Specifically, we find that for appropriate values of the parameters 
of the model, gradient instabilities of vector perturbations can be avoided. 

Outline : The outline of this paper is as follows. In Section 2, we introduce massive bigrav¬ 
ity with a consistent double coupling to matter. In section 3, we assume a flat, homogeneous 
and isotropic universe, and describe some properties of the solutions in two possible different 
branches, which come from two solutions to the Bianchi constraint of the theory. In Section 4, 
we focus on the previously less explored of the two branches. We show specific analytical ap¬ 
proximate solutions for this branch (as introduced in Section 3) in some relevant regimes (early 
times, bouncing period, and late times) for different sets of parameters, as well as numerical solu¬ 
tions. Finally, in Section 5 we summarise our results, discuss consequences and future directions. 

Conventions: Throughout this paper Greek indices such as g, u denote spacetime indices. In¬ 
dices of tensors depending on one metric only will be raised and lowered with their corresponding 
metric, whereas raising and lowering procedures for quantities depending on more than one met¬ 
ric will be explicitly specified where required. The Einstein summation convention is implied 
as usual. In addition, bracketed indices ( g ), and (/), label tensors to clarify their relation 
to/dependence on the two different metrics g /iu and to be considered. These label indices 
are not automatically summed over and whether they are upper or lower indices carries no 
meaning. Finally, we will be using Planck units. 


2 Bigravity model with double matter coupling 

Let us start with the massive bigravity action as proposed in [3]: 

4 

^2 P nen (VS' -1 /) +s m , (2.1) 

71=0 

which includes two Einstein-Hilbert terms for two metrics g and with mass scales given 
by M g and Mf, and their associated Ricci scalars Rt g \ and R(f), respectively. Throughout this 
paper we will use M g = Mf = M p \, where M p \ corresponds to the Planck mass. In addition, 
these two metrics have a very particular interaction term which includes a mass scale m, five 
dimensionless parameters /3 n , and the functions e n (X), which correspond to the elementary 
symmetric polynomials of the eigenvalues of the matrix X = yj g -1 /, which in turns satisfies 
(X 2 )^ = g ,J ' a fau- As previously mentioned, this interaction term is the only way to non- 
dynamically couple these two metrics without introducing a Boulware-Deser ghost. Finally, we 
also include a matter coupling S m of the form 

= f d 4 *^f®£ m [4> 4 , gf u \, (2.2) 


^^2 n n _ n 

S = ~2~ J d * x V-9 R (g) + ^f j d^xy/^fR^-m 4 J d 4 Xy/^g 


where £ m is the matter Lagrangian and all matter fields are minimally coupled to a single 
effective metric c/®® in accordance with the weak equivalence principle. Metric g®® is in general 
a function of both metrics in the theory, g )iu and f pi/ , and it will correspond to the physical 
metric describing the space-time, as from eq. (2.2) we notice that matter will follow geodesics 
of this effective metric. 
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In this paper we focus on the double matter coupling proposed by [70, 71], working in the 
metric formulation used by [70] enabling easier comparison with standard cosmology. In this 
double coupling, the effective metric g®® is given by 3 : 

afv = a2 9fiu + 2af5g /ia X a v + f3 2 f pu , (2.3) 

where a and /3 are two dimensionless arbitrary parameters. Notice that singly-coupled cases 
are recovered by setting either a or /3 to zero. As previously mentioned, massive bigravity with 
matter coupled to the effective metric given by eq. (2.3) does not introduce the Boulware-Deser 
ghost at least up to the strong coupling scale A 3 = (m 3 M p 1 ) 1 / 3 , and thus this theory can safely 
be considered as an effective held theory with a cutoff at A 3 or above. In addition, this double 
coupling was studied in the context of massive gravity (only one dynamical metric), where hat 
FLRW solutions were found [70]. These solutions identified qualitatively new massive gravity 
cosmologies, as no exact hat FLRW solutions are present with single couplings [ 88 ] (although 
solutions with cosmologies arbitrarily close to FLRW are possible). Furthermore, it was found 
that the Boulware-Deser ghost does not propagate at all at the level of linear perturbations 
around FLRW backgrounds [70]. 4 

In the bimetric context, action (2.1) and the double matter coupling with the effective 
metric given by eq. (2.3) have the property of being invariant under the following transformation: 


9gv t-)- j3 n -H- /I 4 -n> Mg •t-t Mj , a 0 /3, (2-4) 

and therefore, contrary to the singly-coupled models (where any potential symmetry between 
the metrics is broken by the matter coupling), both metrics here have a similar role. 

The equations of motion of the bimetric action (2.1) with the previously mentioned double 
matter coupling are [75]: 

V^9G {g)pa (F^g^ + F'V'*) = V Z g^T^a(2^S v a + af^F Xa 5 v p + a/" A F Afr ^), (2.5) 
V^fG {f)pa (F^r + F pv f ap ) = P{2a8 p 5 v a + (3g pX F x j; + f5g uX F x J p p ), (2.6) 

where G^g)^ and G(f) pv are the modified gravitational equations in the absence of matter for 
the metrics g pi , and respectively. Explicitly, G{ g )gu and G(f)g, u are given by: 




9gvR(g)+ 9 £( 1) fin (V 9 V) + 9v (V 9 1 


4 3 




n =0 


4 3 


(2.7) 


~fgu R (f) + — EM) 


n =0 


UFn), ) + UY ( A n)p [FF 1 9 


( 2 . 8 ) 


3 Note that in the vielbein formulation the effective metric vielbein takes a remarkably simple form: a linear 
superposition of the vielbeins for g^ and f^ [71]. 

4 The mentioned double coupling has other interesting features as well. For instance, one-loop corrections from 
matter loops do not detune the particular ghost-free structure of the interaction term of the theory, and thus the 
Boulware-Deser ghost is not re-introduced by these corrections [70, 89]. Note that this result is expected to hold 
for arbitrary matter loop order, as matter loop corrections should sum up to an effective cosmological constant 
as discussed in [71]. 
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where the functions Y^ u (K), written in matrix notation, are given by: 


Y ( o) = 1 , 

y (1) =x - i[X], 

Y(2) =X 2 -X[X] + ^l([X] 2 -[X 2 ]), 

y (3) =X 3 - X 2 [X] + -X ([X ] 2 - [X 2 ]) - -I ([X ] 3 - 3[X] [X 2 ] + 2 [X 3 ]) , (2.9) 

2 6 

where I is the identity matrix and [X] stands for the trace of the matrix X. Here, it is understood 
that (n) is a label taking the values n = 0,1, 2, 3,4. In addition, in eq. (2.5)-(2.6) we have defined 
= 5 /w *X% with inverse F^ v such that F^ a F au = 5„, and the stress-energy tensor T^ u as: 


rp^W _ 2 / 5S m \ 

- V SgfJ ’ 


( 2 . 10 ) 


where it is understood that the indices of T^ are lowered and raised using the effective metric 
g In this sense our T is really a , but we drop the eff label in order to avoid clutter. 

We remark that massive bigravity with double matter coupling is invariant under (a single 
copy of) general coordinate transformations but, contrary to GR, deriving the covariant conser¬ 
vation of the stress-energy tensor with respect to the effective matter metric from the Bianchi 
constraints is not as straightforward. In this paper we will be working with matter (perfect 
fluids, scalar fields) for which this conservation is automatically satisfied, i.e. where 


V|f T***' = 0 


( 2 . 11 ) 


is a direct consequence of the matter equations of motion and where corresponds to the 
covariant derivative with respect to < 7 ®®. For more general types of matter an extra assumption 
may need to be imposed for (2.11) to hold true. The Bianchi constraints will nevertheless 
play a special role, as we will see in the next section. Specifically, in massive bigravity with 
double matter coupling, the Bianchi constraints lead to a consistency equation that identifies 
two qualitatively different branches of solutions. 


3 Cosmological background 

In this section we focus on cosmological predictions of the model presented in the previous sec¬ 
tion. In particular, we show the main equations describing the evolution of a flat homogeneous 
and isotropic universe with a perfect fluid. We also show two possible branches of solutions 
allowed by the Bianchi constraint, as well as some of their properties. 

We assume flat FLRW metrics. Explicitly, both metrics and g^ will be written in 
conformal time r as: 

ds 2 = a 2 [— X 2 dr 2 + S t jdx l dx^] , (3.1) 

ds 2 = a 2 [—X 2 dr 2 + 5ijdx l dx^ . (3.2) 

By means of eq. (2.3), we find that the physical space-time metric < 7 ®® shares the same symme¬ 
tries and therefore can be written in the same form: 

ds 2 s = a 2 (—X 2 dr 2 + 5ijdx l dx , (3.3) 
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where the scale factor a and the shift function X are related to those of the metrics g gi/ and 
fnv by 5 : 


a = aa g + /3a/, 


X = 


aClgXg + fldfXf 

a 


(3.4) 

(3.5) 


Note that the conformal time r will be related to the physical time t by aXdr = dt. From 
eq. (3.4) we find that the comoving Hubble ratio % of the space-time metric g e ^ is related to 
those of the metrics g^ and f gu by: 


aHg + PNU f 

(P N + a) 


(3.6) 


where we have defined N = a // a g , H = a'/a, and Hi = a)/a, (with i = (g. /)). In addition, we 
couple g to a perfect fluid with a stress-energy tensor given by: 


T^u = {p + p)u^u u + pSff, 


(3.7) 


where p is the pressure of the fluid, p its energy density, and u M = ((aX) - 1 , 0 , 0 , 0 ) its isotropic 
4-velocity vector. 

By substituting eq. (3.1)-(3.2) and eq. (3.7) into eq. (2.5)-(2.6) we get the following equa¬ 
tions of motion determining the evolution a flat homogeneous and isotropic universe: 


U 


2 

9 


K 

H'f 


X 2 n al 


9 9 


exp —7 + ?n 4 (/3o + 3N/3i + 3 X 2 /32 + X 5 (3%) 


X' 

2 Ug-f + m 2 g + m A NZa 2 g X g (X f - X g ) + 


a 4 Xg(pXg+pX)(X ~Xf) 


(■ Xf-Xg)a 2 


x / a / 


/3p^- + m 4 (/3iX - 3 + 3X" 2 /3 2 + + ^ 4 ) 


X' 


2H f -^ + 2Hi + 
A " J 


m 


A Za 2 Xf(X g — X/) a i Xf{pXf + pX){X — X g ) 


X 3 


+ 


(X 9 -X/)a 2 


(3.8) 

(3.9) 

(3.10) 

(3.11) 


where we have chosen X = diag(X/X/X 9 , X, X, X) 6 . Note that eq. (3.9) and eq. (3.11) are 
redundant, and can be obtained by taking the derivative of the Friedmann equations (3.8) and 
(3.10), respectively. 

In addition, we will have the matter conservation equation: 


p' = -3 H(p + p), 


(3.12) 


and the following Bianchi constraint: 

{pa 2 a/3 - m A a 2 Z) {X g H f - X f U g ) = 0, (3.13) 

5 The fact that scale factors linearly superpose in this way is a direct consequence of the effective metric 
vielbein being a linear superposition of the vielbeins in the theory. 

6 The massive bigravity interaction term has an ambiguity as the matrix X = \Jg _1 f is not completely 
determined, and a choice must be made. The solution chosen in this paper is the simplest one, and the one that 
can give continuous solutions in the presence of singularities [90, 91]. 
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where we have defined Z = j3\ + 2(3 2 N + foN 2 . As previously mentioned, this last equation plays 
an important role as it divides the possible solutions into two branches. Branch I will be defined 
by the condition (pa 2 a/3 — m 4 a 2 Z) = 0 and Branch II by the condition (X g T~L / — XfT-t g ) = 0. In 
what follows we analyse these two branches, which have qualitatively very different solutions. In 
order to have a clear picture of the physical solutions, and for ease of comparison, this analysis 
will be mainly done using the Friedmann equation for the space-time metric g which can be 
derived from the background equations presented above. In addition, from now on, without loss 
of generality, we choose the time coordinate r such that X = 1. 


3.1 Branch I 

This branch is defined by the following Bianchi constraint: 


(pa 2 a. (3 — m 4 a g Z) = 0. 


(3.14) 


Let us analyse this constraint when p 
rewritten as: 


wp, with w constant. If w ^ 0, 
m 4 Z 


^ (a + (3N) 2 wa/3 


the constraint can be 
(3.15) 


Notice that for w ^ 0, there is a particular choice of parameters that eliminates the N- 
dependence in eq. (3.15): /5/5i = a/3 2 with @ 2 @\ = a 2 @ 3 - However, this choice will be avoided 
as it does not give a valid description of the universe at all times, as it implies p = constant 
always. In any other situation, this equation gives a relation between the ratio N and the energy 
density p of the perfect fluid. 

Using eq. (3.6) and the constraint given by eq. (3.15) we can find the Friedmann equation 
for the space-time metric in this branch in terms of a, N and p. However, expressed in this 
way the Friedmann equation has a rather complicated expression, which we explicitly show in 
Appendix A, eq. (A. 14), together with its derivation. This Friedmann equation shows that, if 
some of the parameters /3 n , a or (3 had opposite signs, in general there would be a divergence 
in % at some finite time during the evolution of the universe, or a violation of positivity of 
one of the Friedmann equations (3.8)-(3.10)'. For this reason, we only consider solutions where 
all parameters have the same sign. Without loss of generality, from now on we assume all the 
parameters to be positive, in which case we will also assume N > 0, as a negative N would 
violate the positivity of 'H 2 , rendering complex the solution for the physical metric. 

Assuming that all parameters are positive, and that a 7 ^ 0 and /3 7 ^ 0, we next analyse the 
possible solutions allowed by the constraint (3.15) in two relevant regimes. 


Early times: We look for possible values of N when p/m 4 00 . As we can see from the 

constraint (3.15), if all parameters are positive, p/m 4 will never reach infinity, but a 
maximum value instead, which will depend on the choice of parameters. This maximum 
value can be found by taking the time derivative of eq. (3.15): 


/ 2m 4 \ [(ap 2 - PPi) + N(ap 3 - pp 2 )] AT , _ , , 

' = -(a + MV 3 - = ( ’ 


(3.16) 


We see that the maximum will be reached when F(N ) = 0 or N' = 0. For the case @3 = 0 
(which will be required later), it is possible to show that N' > 0, except in an asymptotic 

7 Note that there might be very particular choices of parameters that avoid these two problems, but we do not 
explore these cases further as these choices would be very restricted. 
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limit, which actually leads to a minimum of p given by p = 0, as we will corroborate in 
the next section. Then, the maximum of p will be determined by F(N), and it will occur 
when F(Nb ) = 0, where: 


_ («/3 2 - /3/3i) 
/3/3 2 


(3.17) 


where we have set 03 = 0. The maximum of the energy density can be found by evaluating 
the constraint (3.15) in N = Nb'. 


Pmax 


m 4 0| 

wa/3 2 (2a02 — /3/3i) 


(3.18) 


Notice that since N' > 0 and N > 0, early times will in general be described by the period 
where N <C 1, which does not necessarily correspond to the region around Nb , and thus 
around the maximum p max . Therefore, if Nb ~ 1 or Nb 3> 1, we say that the maximum 
Pmax is reached at an intermediate stage, instead of early times. 

The form that the Friedmann equation takes at early times and around Nb will depend 
on the choice of parameters. As such, in the next section, we will show the evolution of 
the scale factor in detail in three different cases: Nb > 0, Nb = 0, and Nb < 0. Finally, we 
notice that since N > 0 we expect to have a maximum of the energy density at N = Nb 
only if Nb > 0. If Nb < 0, the value Nb will never be reached, as the universe would start 
at N = 0 because a negative N would violate the positivity of 'H j. As we will see in the 
next section, in this last case, the maximum value of the energy density will be reached 
when N = 0 and p’{N = 0) / 0. 

Late times: We look for possible values of N when p/m 4 —> 0. From the constraint (3.15) 
we see that the only possibilities are that N —> 0 (if (3± = 0), or N —> 00 (if 03 = 0). 
Notice that these two cases are equivalent due to the transformation given by eq. (2.4). 
Therefore, from now on, without loss of generality, we will assume that 03 = 0, and thus 
N —> 00 when p/m 4 —> 0, which will correspond to late-times. 

Notice that, if in addition we set 0i = 0, we would also have p/m 4 —> 0 when N —> 0. 
Since N' > 0 and N > 0, the limit N —> 0 would correspond to early times, and thus if 
Pi = /?3 = 0 , the universe will bounce, and p/m 4 —> 0 in the infinite past and future. 

If & 7 ^ 0, the Friedmann equation (A. 14) at late times, i.e. for N 1, approximates to: 


n 





(3.19) 


where we recall that we have set X = 1. Equation (3.19) gives a solution that corresponds 
to a de-Sitter phase, with a cosmological constant determined by the parameter 04. 


3.2 Branch II 

Branch II has been previously studied in detail [14, 86, 87], and here we simply summarise some 
relevant results found. This branch is defined by the following Bianchi constraint: 

X g 'Hf-X f 'Hg = 0, (3.20) 
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and when combined with the two Friedmann equations (3.8) and (3.10), we get: 


m 4 [ft + (3/3 2 - 0o )N + 3N 2 (0 3 - 0i) + (0 4 - 30 2 )iV 3 - N 4 f5 3 ] 

+ (0 -aJV)(a + PN) 3 p = 0. (3.21) 


Similarly to Branch I, eq. (3.21) gives a relation between the ratio N and the energy density 
p of the perfect fluid. Using eq. (3.6) and the constraint given by eq. (3.21), we find that the 
Friedmann equation for the space-time metric (the effective matter metric g e g) becomes: 


r 


W = 3 


otpiot + ./V/3) H- 


m 


(a + N/3 ) 2 


(0o + 3-/V01 + 3IV 2 02 + AT 3 /5 3 ) 


(3.22) 


Assuming that a/0 and 0 7 ^ 0, we next analyse the solutions allowed by the constraint 

(3.21) in two relevant regimes. 


Early times: We look for possible values of N when p/m 4 —> 00 . In this case, the constraint 

(3.21) can be satisfied if N —> /3/a or N —> —a/0 8 . As pointed out in [14], the only viable 
situation is when N — > /3/a. In this case the Friedmann equation (3.22) approximates to: 

V 2 = j(a 2 + (3 2 )p, (3.23) 

which corresponds to GR with a modified gravitational constant. 

Late times: We look for possible values of N when p/m 4 —> 0. In this case, the constraint 

(3.21) can be satisfied if AT —>- 0 (if 0i = 0), N —> 00 (if 03 = 0), or N —> N 9 , where 
N is a non-zero constant. Notice that the two first types of solutions are related by the 
transformation given by eq. (2.4), and then they generate the same observables. Therefore, 
without loss of generality, we can only consider the case where N —>• N, where N is now a 
constant including zero. Around the value N = N, we find that the Friedmann equation 

(3.22) approximates to: 

V 2 = —A, (3.24) 

3 

where A is a constant given by: 

4 

A = 7 ™- rm2 (0o + 31V0r + 31V 2 02 + 1 V 3 0 3 ). (3.25) 

(a + N(3) z 


From eq. (3.24) we see that the evolution of a approaches a de-Sitter phase. In the 
particular case of N = 0, the cosmological constant will be given by the parameter 0 q. 

Particular choices of parameters were compared to data for Branch II in [14] , where viable 
(background) solutions for a were described. However, at the level of linear perturbations 
tachyonic, gradient and ghost instabilities were found for tensor, vector and scalar perturbations, 
respectively [ 86 , 87], rendering this branch problematic, as perturbations would either break the 
perturbative approximation, or require fine-tuning. 

s Note that in both cases a particular relation between the parameters a, f) and /3„ has to be assumed in order 

to satisfy the constraint in eq. (3.21). 

9 Note that this particular constant value of N is not arbitrary, but one such that the constraint (3.21) is 
satisfied with p = 0. 
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4 Branch I 


The presence of instabilities at the level of linear perturbations in Branch II motivates searching 
for viable cosmological solutions. In this section we address this problem, at the level of the 
background, by analysing the evolution of the space-time metric (3.3) in Branch I in detail. 
In particular, we distinguish three different choices of parameters, according to their values of 
IV;,, and find approximate analytical solutions in some relevant regimes, as well as numerical 
solutions. As previously mentioned, we will be using the Friedmann equation for the physical 
metric given by eq. (A. 14) for analysing the behaviour of the scale factor a in this branch, instead 
of the Friedmann equations for each metric and f^ u . We do this for ease of comparison 
with other models, and as it is physically clear. As explained in Appendix A, eq. (A. 14) can 
be written in terms of N only, and therefore its behaviour can be simply analysed in terms of 
N, which is what we do throughout this section. We recall that we will be assuming fa = 0, in 
order to have p/m 4 —> 0 at late-times. 

4.1 Positive case: Nj, > 0 

In order to have IV& > 0, we could have (5\ = 0, or 0\ / 0 and 02 / 0 such that {a.02 — 00i) > 0. 
In what follows, we find the evolution of the scale factor a during the radiation-dominated era, 
in three relevant regimes: early times (N <C 1), intermediate bouncing regime (N ~ A),), and 
late times (N 3> 1). 


Assuming w = 1/3, during early times, i.e. for N <C 1, eq. (A. 14) approximates to: 


T~L = — Ha ; H = m 


2 (a/32 ~ PPi) V00o + 3aPi 
(a/32 + /3/3i) a^30 


> 0 , 


(4.1) 


and if 0 \ and 0 o are not zero at the same time, the evolution of the scale factor in 
conformal time will be given by a = 1 /(Hr), which translates into the following solution 
in the physical time t: 

a(t) = ae~ Ht . (4.2) 

Therefore, the universe is contracting, yet still accelerating with a > 0. Here, the value 
of a is determined by initial conditions. From eq. (3.15) we can also find the evolution on 
N at early times: 

P ~ [«/3i + 2(a/3 2 - PPi)N] , (4.3) 

and then N as & function of physical time is given by: 


m = 


a(po r a 3 0(ma) 4 e 4Ht — 30i) 

6(a02 ~ 00i) 


(4.4) 


where we have used that p = po r /a 4 , where po r is the energy-density of radiation today. 
As we can see from this equation, if 0\ / 0, as t decreases, N decreases and the value 
N = 0 is reached at a finite time to, which gives a non-zero value of the scale factor 
a{t 0 ) = ae~ Ht °. The universe ends at this time, as otherwise N would take negative 
values and thus Ti'j would become negative, rendering the solutions complex. 

On the other hand, if 0\ = 0, according to eq. (4.4), the value N = 0 is reached in the 
infinite past t —> —oo, limit at which a —> oo. Thus, this case has an infinite contracting 
past, and no violation of positivity is present. For this reason, from now on we focus on 
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the /? 2 -only case. We remark that when t —> —oo, N = 0 and N' = 0. This asymptotic 
limit is the only one that gives p' = 0 and N' = 0. 

Finally, we analyse the next-to leading order term of eq. (4.1). When j3\ = 0, the Fried¬ 
mann equation approximates to: 


n 


—Ha + CNa] C 


1 \/3m 2 (A/3 2 — 2>f3^ 2 a 2 + (3 2 (3^ 2 ) 

3 yffa A)a 2 /3 


(4.5) 


Since from eq. (4.3) we see that p oc IV if we choose parameters such that C < 0, the 
evolution of the universe would mimic GR with a cosmological constant (determined by 
/3o) at early times, but in a contracting universe. If we set /3p = 0, the universe would 
approach GR without a cosmological constant (in a contracting universe, as before), as 
the Friedmann equation would approximate to: 


% = —m 2 



a = —a 



(4.6) 


where in the last step we used eq. (4.3). 

• In what follows, we analyse the evolution of the scale factor in the intermediate regime 
given by N ~ Nf, = a//3. As previously explained, at this time the universe reaches a 
maximum of p, and therefore it bounces. Assuming w = 1/3, for N ~ 7V&, the Friedmann 
eq. (A. 14) approximates to: 


H = Ba(N-N b y, 


HgbHfgP^ 

4a(/3H g b + otHf b ) ’ 


(4.7) 


where we have defined H g and Hf such that X 2 a 2 H 2 = H 2 and X'ja 2 Hj = H 2 , and the 
subindex b denotes evaluation at the bouncing time. In order to find the solution for a, we 
rewrite eq. (4.7) in terms of a only. From the constraint (3.15) we find that the relation 
between a and N around the bounce is given by: 


a = o 0 + ai(N - N b ) 2 ; a 0 


2por-a 2 /3 2 l 1/4 _ ap/3 2 

3/?2^ 4 ’ 1 16a 2 ’ 


(4.8) 


where por is the energy-density of radiation today. Combining this last equation with 
eq. (4.7), the Friedmann equation becomes: 


H = Ba 0 J (4.9) 

V ai 

whose solution in physical time is given by: 

a(t) = a 0 + - t 0 ) 2 . (4.10) 

4ai 

As expected, there is a quadratic bounce in this regime, at which the scale factor takes 
its minimum value a(t = to) = ap given by eq. (4.8), which is completely determined by 
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the parameters. At this minimum a o, the energy density reaches a maximum given by 
eq. (3.18) with /3i = 0: 


3 m 4 / 

Pmax= 2 ( 02 ^ 2 ' 

Around the bounce, from eq. (4.8) we find that N evolves as: 


(4.11) 


N(t)=N b + -±(t-t 0 ). (4.12) 

2 ai 

Finally, we remark that if j3\ ^ 0, the overall evolution for a and N during the bouncing 
period would have been the same. 


• Next, we analyse the evolution during late times, i.e. for N 3> 1. In eq. (3.19) we already 
found that the universe approached a de-Sitter phase determined by /3 4 at late times, and 
then 


a(t) = ae Hot , 


(4.13) 


where a is determined by initial conditions and Hq is given in eq. (3.19). 
constraint (3.15) we find that: 

6 m 4 /?2 
^ a(3 3 N ’ 

which means that N evolves as 


From the 
(4.14) 


m 


6m 4 /32 

a/3 3 p 0r 


a 4 e 4Hat , 


(4.15) 


where po r is the energy-density of radiation today. 

Next, we analyse the next-to leading order term in the large N approximation: 


T~i = m 2 


I Pa C 

W a+ N a] 


C = - 


m 2 [(4/l 4 a 2 - 3/3 2 /3 2 )\/^2 + a 2 /?^ 2 ] 
a(3 2 y/ 3/3 4 /?2 


(4.16) 


But by using eq. (4.14) the Friedmann equation can be rewritten 


% = m z 


I Pi CaP 3 

W a+ ^% pa - 


(4.17) 


From here we see that if C > 0, at late times the evolution approaches GR with a 
modified gravitational constant and a cosmological constant. Similarly to the evolution 
at very early times, if /3 4 = 0, the solution mimics GR without a cosmological constant, 
as the Friedmann equation approximates to: 


n = = ^ <4 ' 18) 

For this reason, we will require Pa ^ 0 to obtain a cosmological evolution with late-time 
acceleration. Finally, we remark that Pi ^ 0, the overall evolution for a and N during 
late times would have been the same. 


We now discuss numerical solutions for a, N, a g , a,f during the radiation-dominated era, 
as shown in in Figure 1. On the left-hand side of that figure, we show the scale factor a as a 
function of physical time t, while on the right-hand side we show the ratio of scale factors N. 
We use arbitrary initial conditions and parameters such that N b > 0. 
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Figure 1 . Scale factor a (LHS) and ratio N (RHS) as a function of physical time t during the radiation- 
dominated era for a = 8 = 1, ,3-2 = 16, and /?4 = 10, flo = 9, m = 1. In this case Nb > 0 and the bounce 
ocuurs when N = 1. 


As in the analytical solutions, in Figure 1 we see that N is a growing monotonic function, 
where N <C 1 characterises early times, and N 1 characterises late times. In addition, 
we observe the intermediate bouncing regime, with a minimum for the scale factor. In this 
numerical example we set /3o 7 ^ 0 , which generates a contracting accelerated period for a during 
early times, and /I 4 7 ^ 0, which generates an expanding accelerated period at late times. In 
Fig. 2 we show numerical solutions for the evolution of scale factors a g and a 7 as a function of 
physical time t during the radiation-dominated era, for the same initial conditions and choice 
of parameters as in Fig. 1. For this bouncing solution, we see that the scale factor a g always 
decays in time, whereas 0 / always grows in time. In fact, this behaviour can be derived from 
the previous analytical solutions found for a and N. We find that at early times a g oc e~ Ht and 
cif oc e 3Ht (when ft 7 ^ 0 ), and therefore a g —» 00 and aj —> (1 in the infinite asymptotic past. 
During late times, a g oc e~ 3Hot and a/ oc e Hot (when ^ 7 ^ 0), and thus a g —> 0 and aj —> 00 in 
the infinite asymptotic future. 

In summary, when Nj, > 0, the f5\ and @2 case differs from the /? 2 -only case only during 
early times, as the former has a finite past with a non-zero value of the scale factor, whereas 
the latter has an infinite past with a —> 00 for t —> — 00 . In the / 32 -only case we found that the 
universe mimics GR with a cosmological constant at early times if /3q 7 ^ 0. In both cases there 
will be an intermediate bouncing regime with a non-zero minimum value of the scale factor, and 
a GR-like regime with cosmological constant at late times if (3^ 7 ^ 0. We also remark that, even 
though some of the solutions found did violate the positivity of 'Hj? at some finite time, they 
could still be viable, as that time would represent the start of the evolution of the universe. 
This is different to what we would find if some of the parameters f3 n , a and /3 had different sign, 
where violations of positivity would be likely to happen at some time in the future, and hence 
the universe would end at a finite time in future. As previously mentioned, we discard this last 
kind of solutions as we restrain ourselves to only look for standard evolutions with an infinite 
future here. 
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t t 

Figure 2. Scale factor a g (LHS) and a/ (RHS) as a function of physical time t during the radiation- 
dominated era for a = /3 = 1, = 20, and Pi = 10, Po = 9, m = 1. In this case Nb > 0. 


4.2 Null case: Nb = 0 

In this case Pi and (3 2 need to be different from zero, and chosen such that ( ap 2 — PPi) = 0. 
As previously explained, we expect a minimum of the scale factor a (and a maximum of p ) to 
be reached when N = 0. Notice that, in this case, early times iV C 1 coincide with the period 
N ^ Nb = 0. Thus, there is no intermediate regime, and in what follows we analyse only two 
relevant regimes during the radiation-dominated era, namely early times and late times. 


At early times, i.e. JV<1, the Friedmann equation (A. 14) approximates to: 


H = rut 


2 //3(3/3i a + fifio) 


12 a 4 


Na. 


(4.19) 


In order to find the solution for a, we rewrite this equation in terms of a only. From the 
constraint (3.15) we find that: 


P 


3m 4 /3i 

a 3 /3 


1 - 

cr 


and thus in this regime the explicit relation between a and N is given by: 

(3 2 a 0 


a = oq + a\N ] ao = 


„ ^,2 a 2 "| 1/4 
P0r(% P 


3/^2 m 4 


ai = 


4a 2 ’ 


(4.20) 


(4.21) 


where po r is the energy-density of radiation today. By means of eq. (4.21), eq. (4.19) can 
be rewritten as: 


H = m 


2 P{3PiOt + PPo) 


12a A a\ 

whose solution in physical time is given by: 

a 0 m 4 (3p±a + PPq) 


ao\/a — ao, 


a(t) = a 0 + 




12a 2 




(4.22) 


(4.23) 
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where to is determined by initial conditions. As expected, the early time period corre¬ 
sponds to the minimum a period. Here, the scale factor evolves quadratically and reaches 
a minimum value ao, given in eq. (4.21). 

In addition, from eq. (4.21), we find the evolution for N as a function of the physical time: 


N(t) = 


?n 2 v //3(3/3ia + f3/3 0 ) 
V3/3 2 


(t -1 0 ), 


(4.24) 


which means that the value N = 0 is reached at a finite time to, which is the moment at 
which the evolution starts as N cannot become negative (as before, this would violate the 
positivity of 77 2 ). In this case, the maximum value of the energy density will be given by 
eq. (3.18): 

_ 3m4 / 3 2 tA 

Pmax ~ a 2p2 ■ (4.25) 


• At late times, i.e. N 1, the solutions of a and N are the same as for the previous case 
with Nb > 0, given by eq. (4.17) if /?4 / 0, and by eq. (4.18) if = 0. Therefore, the 
late-time universe will have a cosmological constant when /J 4 7 ^ 0 . 


Next, in Figure 3, we show numerical solutions during the radiation-dominated era and 
when Nb = 0. On the left-hand side we show the scale factor as a function of physical time 
t, while on the right-hand side we show the ratio of scale factors N. We use arbitrary initial 
conditions and a choice of parameters such that Nb = 0. As seen in the analytical solutions, N is 
a growing monotonic function, where N -C 1 characterises early times and N 1 characterises 
late times. In addition, we have a finite past of the universe, where the scale factor a starts at 
t = —2.5 at a non-zero minimum value. Since we set [3^ 7 ^ 0 we observe a late-time accelerated 
expansion of the universe. 




t t 

Figure 3. Scale factor a (LHS) and ratio N (RHS) as a function of physical time t during the 
radiation-dominated era for a = (3 = 1, = Pi = 1, and (3 4 = 10, /3q = 0, m = 1. In this case Nb = 0. 


Furthermore, Figure 4 shows the evolution of the scale factor a g and aj as a function of the 
physical time t during the radiation-era. We chose the same initial conditions and parameter 
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Figure 4. Scale factor a g (LHS) and a/ (RHS) as a function of the physical time during radiation- 
dominated era for a = /3 = 1, = Pi = 1, and /J 4 = 10, (3q = 0, m = 1. In this case Nb = 0. 


values as in Fig. 3. In Fig. 4 we see that a g always decays in time, whereas a/ always grows, 
similarly to the Nb > 0 case. This behaviour can be derived from the previous analytical 
solutions found for a and N. We find that during early times a g — a gmax oc — (t — to) and 
aj oc (t — to), where a 5max is the maximum value of a g at the start of the evolution, at t = to- 
Then, at the beginning of the evolution a g = a gmSLX and aj = 0. On the contrary, during late- 
tirnes a g oc e^ 3Hot and a/ oc e Hot (if /?4 / 0), and therefore a g —> 0 and a,f —> 00 in the infinite 
future. 

Summarising, when Nb = 0, we need ft / 0 and fa / 0 such that (a/3 2 — /3/3i) = 0. 
The universe starts at a finite time to at which there is a minimum non-zero value of the scale 
factor ao, completely determined by the parameters (see eq. (4.21)). The universe expands 
quadratically at early times, and approaches GR with a cosmological constant at late times 
when /?4 / 0. 

4.3 Negative case: Nb < 0 

In this case we could have (3 2 = 0, or /?2 / 0 and f3\ 7 ^ (J such that (a02 ~ /ftSi) < 0. As 
previously mentioned, there will be no bounce, as N cannot be negative, and thus Nb will never 
be reached. Thus, similarly to the previous case with Nb = 0, we will analyse only two relevant 
regimes during the radiation-dominated era, namely early times and late times. 

• During early times, i.e. JV<Cl, the Friedmann equation is the same as eq. (4.1), but now 
H < 0. We can then write the Friedmann equation as: 

U = \H\a , (4.26) 

and therefore the evolution of the scale factor in physical time is given by: 

a(t) = a 0 e mt - to \ (4.27) 

where to is determined by initial conditions. We can see that during early times there is an 
exponential growth (regardless of the value of /3q). In what follows we find the evolution 
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for N also. From the constraint (3.15) we find the relation between a and N to be: 


CL — fflQ T fli-A; ®o 


Por« 3 /3 

3/Sim 4 



a i — ao 


{(5Pi ~ a fa) 

2a(3\ 


> 0 , 


(4.28) 


where po r is the energy-density of radiation today. Then, the evolution of N in physical 
time is the given by: 

N(t) = ^ (e |J?l(t_to) - l) , (4.29) 

where we can see that N grows exponentially in time, and N = 0 is reached at a finite 
time to, when the scale factor a reaches its minimum value ao- At this time to the universe 
starts as N cannot be negative (otherwise it would violate the positivity of Ti'j as before). 
The maximum value for the energy density will given by eq. (3.15) evaluated at N = 0: 


3m 4 /3i 

Pmax = 7“T“ 

a 6 p 


(4.30) 


Notice that, contrary to the previous cases, at t = to, p is at its maximum, but p'(to ) 7 ^ 0. 
Finally, we mention that the next-to-leading order term in eq. (4.26) is proportional to 
y/N. Since N is related to p at early times through eq. (4.3), this correction term goes as 
y/ Pmax P■ 


• At late times, i.e. N 1, the Friedmann equation (A. 14) approximates to eq. (4.17) when 
/?2 / 0 and /?4 7 ^ 0, and eq. (4.18) when /?2 / 0 and /?4 = 0. On the other hand, if fa = 0, 
eq. (A. 14) approximates to: 


n 



1 m? \/2oiP^ a 
6 fayffa 7F 


(4.31) 


However, from the constraint (3.15) we find that: 


2>m A Pi 1 

a/3 3 W' 


(4.32) 


and then the next-to-leading order term in eq. (4.31) evolves as p 1 / 4 . Therefore, if fa = 0, 
a does not mimic GR, even when ^4 7 ^ 0. If in addition fa = 0, then 


H = m 2 




(4.33) 


and the solution asymptotically approaches GR without a cosmological constant at late 
times. 


A numerical solution during the radiation-dominated era when Nf, < 0 is shown in Figure 
5. On the left-hand side we show the scale factor a as a function of physical time t , while on 
the right-hand side we show the ratio of scale factors N. We use arbitrary initial conditions 
and parameters such that Nj, < 0. 

As seen in the analytical solutions, in Fig 5 we see that N is a growing monotonic function, 
where JV <C 1 characterises early times, and N 1 characterises late times. In this solution 
the universe starts with a non-zero minimum value of the scale factor, with an accelerated 
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Figure 5. Scale factor a (LHS) and ratio N (RHS) as a function of the physical time during radiation- 
dominated era for a = P = 1, P 2 = 1, Pi = 10, and P 4 = 10, /3q = 10, m = 1. In this case < 0. 


early period. In addition, since we set /?4 ^ 0, we also have a late-time accelerated period. 
Furthermore, Figure 6 shows the evolution of the scale factor a g and a/ as a function of the 
physical time during the radiation-dominated era. We set the same initial conditions and 
parameter values as in Fig. 5. In Fig. 6 we see that a g always decays in time, whereas aj always 
grows, similarly to the previous cases. In fact, this behaviour can be derived from the previous 
analytical solutions found for a and N. We find that during early times the evolution is such 
that (a g — a gmax ) oc (1 — and a,f oc (el-^K* - * 0 ) — 1), where at to the universe starts. 

Thus, a g = a g max and a/ = 0 at the beginning of the evolution. During late-times a g ~ e ~ 3 H o * 
and af ~ e Hot (if ,64 7^ 0 and ( 5 2 ^ 0 ), and therefore a g — > 0 and a/ —> 00 in the infinite future. 



t t 


Figure 6. Scale factor a g (LHS) and a/ (RHS) as a function of the physical time during radiation- 
dominated era for a = ft = 1 , /3 2 = 1 , Pi = 10, and P 4 = 10, Po = 10, m = 1 . In this case Nj, < 0. 
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Summarising, when iV& < 0, the fa and fa case differs from the /?i-only case during late 
times only, as the former can mimic GR at those times, if fa ^ 0, whereas the latter cannot. 
During early times the evolution in both cases is the same, namely a finite past with a minimum 
scale factor ao, and an exponential growth of a in physical time. 

4.4 Radiation and dust 

Previously, we studied the evolution of the universe in Branch I during the radiation-dominated 
era only. In what follows we include the effect of dust, focusing on the bouncing fa-oxAy solution. 
We emphasise that we do not assume radiation to be completely negligible as it is usually done. 
This is because radiation plays a special role through the constraint given by eq. (3.14), as this 
constraint depends solely on the pressure of the perfect fluid, hence radiation should not be 
fully neglected (see Appendix B). When adding dust, the Friedmann equation will be given by 
eq. (A. 14), where the explicit dependence on p and w will correspond to the values for radiation, 
while the implicit dependence on p (via the functions H g and Hf defined in the appendix), will 
correspond to the value for dust and radiation. This is due to the fact that in the derivation of 
eq. (A. 14) all its explicit dependence on matter came from the constraint (3.14), and therefore 
it only includes fluids with non-zero pressure. After making these replacements, we analyse the 
Friedmann equation in three relevant regimes, namely early times, bouncing period, and late 
times. 


• At early times, i.e. when N <C 1, the Friedmann equation (A.14) approximates to: 


H = - ro2 vH“ ■ 2^ m Pma • (434) 

where p m is the energy-density of dust. As we can see, the leading order term is the same 
as the one found during the radiation-dominated era in eq. (4.1) with fa = 0. However, 
the next-to-leading order term differs, but the solution still asymptotically approaches GR 
(in a contracting universe) at early times. Similarly to the solution during the radiation- 
dominated era, if we set fa = 0, the Friedmann equation would be of the form: 


T~L = —a 



(4.35) 


and the solution would approach GR without a cosmological constant at early times. 


• Around the bounce, i.e. when N ~ Nb, the Friedmann equation is given by eq. (4.7), 
where the only difference is that previously the explicit expressions for H g b and Hjb in 
the coefficient B included radiation only, but now they include radiation and dust. In 
addition, relation (4.8) also holds when dust is included, and thus the expression for the 
minimum value of the scale factor does not change, and so it will depend on the radiation 
energy-density only. 

• During late times, i.e. when N 1, the Friedmann equation (A. 14) approximates to: 

M = m2 /^° + 2^73R A ” a ’ <436> 

where p m is the dust energy-density. Thus, the solution approximates GR with a cos¬ 
mological constant at late times when fa ^ 0. Similarly to the results found during the 
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radiation-dominated era, if /h = 0 the solution does not have a cosmological constant but 
still approaches GR (without a cosmological constant) at late times, as the Friedmann 
equation approximates to: 


U 



(4.37) 


Next, we show numerical solutions when radiation and dust are included. Figure 7 on the 
left-hand side shows the evolution of the scale factor a as a function of physical time t, while 
on the right-hide side it shows the evolution of the ratio of scale factors N. We use arbitrary 
initial conditions and a choice of parameters such that > 0 . We observe the same overall 
behaviour as in Fig. 2. Here we have chosen the parameters in a way that the evolution of the 
scale factor a presents prolonged decelerating periods just before and after the bounce, which 
generate a difference between this figure and Fig. 2. During these periods, the cosmological 
constants /3q and /J 4 do not yet dominate. 




Figure 7. Scale factor a (LHS) and N (RHS) as a function of the physical time t during radiation 
and matter dominated eras, when /3 = a = 10 -1 , /?4 = 3 x 10 -10 , /3q = 2 x 10~ 10 , fa = 5 x 10 4 , and 
/?3 = /3i = 0. In this case 1V& > 0. 


5 Summary, discussion and conclusions 

In this paper we analysed a new branch of flat homogeneous and isotropic solutions in massive 
bigravity with double matter coupling. In this branch, even though generic choices of parameters 
can lead to divergences of the Hubble rate T~L at finite times and/or violations of positivity of the 
equations, the family of parameters j3 n , a, (5 > 0 gives viable cosmological background solutions. 
All these solutions share the following characteristics: they have a minimum non-zero value of 
the scale factor, an accelerated period during the radiation-dominated era, and an accelerated 
period at late times if /?4 7 ^ 0. Furthermore, all the solutions have a finite past, except when 
Pi = @3 = 0. This last case is the only one that avoids divergences and violations of positivity 
at all times, and thus gives a universe with an infinite past and future. In Table 1 we summarise 
the main solutions found in this branch. 
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Non-zero di, 2,3 

Past 

Intermediate 

Future 

N b > 0 

P1-P2 

P2 

Finite exponential/polynomial contraction 
Infinite exponential/polynomial contraction 

Bounce 

Bounce 

GR+A 

GR+A 

O 

II 

£ 

P1-P2 

Finite quadratic expansion 

- 

GR+A 

N b < 0 

P1-P2 

Pi 

Finite exponential expansion 
Finite exponential expansion 


GR+A 
not GR+A 


Table 1. Summary of evolutions of the scale factor a in Branch I when /3 n ,a,f3 > 0. In the first 
column, solutions were separated in three categories according to the value of N b (see eq. (3.17)), namely 
when Nt > 0, N b = 0, and N b < 0. We recall that TV), corresponds to the value of the ratio of scale 
factors N = af/a g at the point where the energy density satisfies p' = 0 (and reaches its maximum 
Pmax for N b > 0, or Nb = 0). In the second column we mention the non-zero interaction parameters 
(/3i, /3 2 , and ^ 3 ). Notice that it is always assumed that a / 0 and /3 7 ^ 0. In the three cases we set 
p 3 = 0 as this gives an infinitely expanding universe in the future. In the rest of the table we describe the 
evolution of the universe in the past and future, and during an intermediate bouncing stage, if present. 
All solutions have a finite past except when Nb > 0 and 7 ^ 0, where there is an infinite contracting 
past. In addition, during the past, the contraction can be exponential (j3o 7 ^ 0 or /3i 7 ^ 0) or polynomial 
(do = /3i = 0) as a function of the physical time t when Nb > 0. For N f) = 0, there will always be a 
quadratic expansion, whereas for Nb < 0 there will always be an exponential expansion. On the other 
hand, when Nb > 0, solutions present a bouncing intermediate stage. Finally, at late times all solutions 
asymptotically approach GR with a cosmological constant if /I 4 7 ^ 0, except when Nb < 0 and /3i 7 ^ 0, 
where the solutions do not mimic GR but do have a cosmological constant. 


The solutions we find are interesting in both early and late universe settings. In an early 
universe context, they naturally lead to an early-time inflating period, without the assumption 
of a fundamental yet unknown field with very constrained self-interactions driving this early- 
time accelerated expansion. In addition, since the minimum value of the scale factor is not zero 
(and the energy density does not diverge), these types of solutions also remove the physical 
Big Bang singularity present in GR. All these features are new in massive bigravity with both 
single and double matter coupling, as all previously found viable cosmological solutions started 
at a = 0 in an early decelerating epoch filled with radiation. However, it is important to notice 
that massive bigravity with double matter coupling is an effective field theory with an unknown 
cutoff, at least equal or higher than its strong coupling scale given by A 3 . Therefore, for the 
solutions found to be viable alternatives to standard inflation, the parameters should be chosen 
in a way that the maximum p max happens at a lower energy scale than the cutoff scale of the 
theory, otherwise the inflating and/or bouncing feature cannot be trusted. For the usual choice 
of parameters, i.e. rri ~ 10~ 33 eV, the scale A 3 would be of the order of 10“ 13 eV, and thus if this 
was the cutoff scale of the theory, no early time evolution found could be trusted. However, in 
[92] an upper bound for the cutoff of the theory was found, given by [(m 3 M^) / (ot/3)] 1 / 5 . If this 
was the actual cutoff, for instance for the bouncing solution with /?i = /I 3 = 0, if rn r^j 10 33 eV 
the bounce could in principle be within the regime of validity of the theory if it happened at 
an energy scale of around 10 7 eV 10 . As this energy scale is approximately that of Big Bang 
nucleosynthesis, we would like the cutoff and the bouncing scale to be higher, which could be 

ll) Here we have considered a/3 ~ 10 -80 and +2 ~ 1 in order to get pj/ax ~ 10 7 eV, according to eq. (4.11). For 
larger values of a/3 the cutoff and bouncing scales would be smaller than 10'eV, and thus BBN could not take 
place. For smaller values of a/3 the cutoff and bouncing scales would be larger than 10 7 eV, but the bouncing 
scale would be larger than the cutoff scale, so the bouncing feature would be outside the regime of validity of 
the effective held theory. This choice of a, (3 is of course fine-tuned - an examination of the naturalness of the 
parameter values chosen here is beyond the scope of this paper, however. 
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achieved by considering a larger m and a larger a/3. However, the tuning on the parameter 
m (and ct, (3) must be done carefully, as it cannot be arbitrarily large, if massive bigravity is 
to satisfy solar system gravity constraints as well. Clearly, a detailed analysis of observational 
and technical naturalness constraints is required to know if (and in what sense) solutions in 
Branch I could be alternatives to standard inflation. Regardless of whether this is the case 
or not, we emphasise that these solutions are of relevance in a late-universe context. As we 
have shown throughout this paper, they provide consistent cosmological evolutions at those late 
times, closely mimicking GR with a cosmological constant in most cases. 

Further work is necessary to determine the full viability of Branch I solutions beyond the 
background behaviour studied in this paper. Here we have analysed their main properties, but a 
detailed study contrasting the solutions we found with data should be carried out. Furthermore, 
it is necessary to analyse the stability and evolution of cosmological perturbations around these 
backgrounds. A general analysis of linear perturbations in Branch I and II was previously 
performed in [87]. Following the standard SVT decomposition [93] and the results of [87], in 
[94] we argue that it is possible to avoid ghost and gradient instabilities on tensor and vector 
perturbations in Branch I, but tachyon instabilities will be present in both. In particular, for 
the bouncing solution with 0i = 03 = 0, ghost and gradient instabilities can be avoided when 
/?4 > (32, 0o > P 2 and 04a 2 + 150 2 02 > 0o0 2 + 1502a 2 . Due to the involved expressions describing 
scalar perturbations, we do not comment on their stability further here. Instead, we leave this 
issue for future work. 

We finally remark that, within the bouncing solutions with 0i = (3% = 0, we might find 
partially massless bigravity theories [95, 96], defined by the following choice of parameters: 

0o = 30 2 = 04, 01 = 03 = o. (5.1) 

Arguably the most interesting property of these theories is that around a de-Sitter space (both 
metrics are de-Sitter 11 ) they display an extra gauge symmetry (besides diffeomorphism invari¬ 
ance) that removes the helicity-0 mode of the massive graviton at linear order (and beyond, 
should it turn out to be a full gauge symmetry). Due to the extra symmetry, these theories prop¬ 
agate 6 DoF, as opposed to the 7 DoF propagated in massive bigravity. If this gauge symmetry 
can be non-linearly completed, these models would provide a very interesting modification to 
GR. It is not known if this gauge symmetry appears around the backgrounds found in Branch 
I, but this should become clear when studying linear perturbations around such backgrounds 
in the future. 
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A Branch I: Friedmann equation for effective metric 

In this section we find an expression for the Friedmann equation for the scale factor a of the 
physical metric g in terms of the functions a, N and p. We recall that we will be setting 
X = 1. 

11 This is different to the late-time de-Sitter phase found in Branch I, where the effective metric expanded 
exponentially in time, but the individual metrics / M „ and did not. In fact, we recall that at late-times the 
individual scale factors evolved as a g oc e~ JHot and a/ oc e Hot . 
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First, assuming w 7 ^ 0, we define the intermediate function F(N ) such that: 


p' = F(N)N', 


(A.l) 


which can be obtained by taking the derivative of the constraint (3.15) defining Branch I. The 
explicit expression for F(N) is given by: 



[(af3 2 - (3f3 1 ) + N(a/3 3 - /3/3 2 )] 
(a + /31V) 3 


(A.2) 


On the other hand, by definition, we also have that: 


N' = N(H f -H g ). 


(A.3) 


Combining this last equation with eq. (A.l) and the conservation eq. (3.12) we get: 


H f = H g 


371(1 + w)p 

FNTig 


(A.4) 


where we have also used that p = wp. By replacing this last equation into eq. (3.6) we find: 


V = V (« + j3N)F 

9 (a + /3N)F + 3f3p(l + w) 


(A.5) 


Next, we would like to eliminate the explicit dependence on 7i g from the last equation, in favor 
of a, N and p. In order to do this we define H 1; for i = (g, /) as: 


Hi = \JHf/(ciiXi) 2 > 0, (A. 6 ) 

and thus 

Hi = ± HidiXi , (A.7) 

where there are in principle four combinations of signs to consider, but only two are relevant: 
if both Tii have the same sign, or different signs. Without loss of generality we consider only: 

Tig = ± HgdgXg, Hf = HfCLfXf, (A. 8 ) 

and the other two cases will be different by an overall sign in 7~L only. Explicitly, from the 
Friedmann equations for g^ and given by eq. (3.8) and eq. (3.10), we get: 


Hg = \l 3 [ap(a + /3A0 3 + m 4 (/3 0 + 3iVft + 3Ai 2 /3 2 + !V 3 /3 3 )], 


H f = \U 


?P ^r + m A (h A" 3 + 3N~ 2 h + 3N-ifo + &) 


(A.9) 


When replacing eq. (A. 8 ) into eq. (A.5), we get two possible expressions for 7i, namely 77+ and 
7i-. 


7l± = ±X, 


H g F a 


(a + j3N)F + 3(5p(l + w) 


Finally, we express X g in terms of a, N and p. From eq. (3.5) we have that: 


* 9 = 


(a + j3N) - (3NX f 


a 


(A.10) 


(A.ll) 
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but, eq. (A.4) can be rewritten as 


n f _ N ( X f \ - lT 3U±(l + w)p 
Ug H g \XgJ FNH g a(a + PN)X g 


(A.12) 


where we have used eq. (A.8) and eq. (3.4). From here we can work out Xj as a function of X g 
and the other background functions. When this result is replaced back into eq. (A.11) we get 
for X g \ 


x 9 = 


(a + j3N)FHfN + 3U±Pp(l + w)(a + f3N)a 


-l 


FN 


caHf A /3Hg 


(A.13) 


Finally, combining this last equation with eq. (A. 10), we find that the final Friedmann equation 
for a is given by: 


n± 


H g H f aFN(a + /3N ) 

±aNH f [(a + (3N)F + 3/3p(l + «;)] + (3H g [N{a + f3N)F - 3ap(l + w)] ’ 


(A.14) 


which corresponds to the Friedmann equation for the physical metric c/®® as a function of p, a 
and N (provided eq. (A.2) and eq. (A. 9)). 

From eq. (A.14) we notice that the dependence of 7i± on p can be removed entirely by 
means of the constraint given by eq. (3.15), in which case we would get: 


H± = aH±, 


(A.15) 


where H± would be a function of N only, and it would represent the Hubble rate in physical 
time. This last way of expressing H± will be particularly useful throughout the paper, as a 
solution for a can be found by analysing the behaviour of N only. 

It is useful to show the explicit form of the brackets in the denominator as a function of 
N only: 

4 

777 

(a + /3N)F + 3/3p(l + w) = a ^ w ^ a [(^ a ^2 + /5/5i(1 + 3 w)) 

+2N {afc + A5 2 (2 + 3w)) + 3/3(1 + w)foN 2 ] , (A.16) 

4 

777 / 

N(a + /3N)F - 3ap(l + w) = . - n9 [-3a/3i(l + w) - 2N (/3/3i + a/3 2 (2 + 3 w)) 

apwya + pJS) z 

—N 2 (2/3/3 2 + a^ 3 (l + 3w))] , (A.17) 

where we have used the constraint given by eq. (3.15). From these last two expressions and 
eq. (A.14) we can see that, if some of the parameters /3 n , a or /3 were negative and some positive, 
it would be likely to have a zero in the denominator for a finite value of N, which in general 
would generate a divergence in H at a finite time, making the solution for a unviable. Even 
though there might be specific situations in which the parameters are chosen such that the 
divergence is avoided (with a compensating zero in the numerator), or the divergence occurs in 
the infinite future, we do not consider those highly specific cases further here. In addition, for 
negative values of the parameters there might also be violations of positivity of H g and Hf at 
finite times, which would render those solutions unphysical as the Friedmann equations for the 
metrics g^ and f gv would become complex. However, if instead all the parameters /3 n , a and 
/3 have the same sign, we can avoid a zero in the denominator, if w > 2/3 and if we chose 
This is the situation we consider in this paper. In particular, and without loss of generality, we 
assume all the parameters a, f3 and /3 n to be positive. 
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B Matter-dominated universe 


In a matter-dominated universe the constraint given by eq. (3.14) would become: 

a? g Z = 0, (B.l) 

where we have set p = 0. This constraint can only be satisfied if a y = 0 (N —> oo), or Z = 0. 

However, since we are assuming all parameters to be positive, Z = 0 can only be satisfied if 

N = 0 when j3\ = 0. According to the previous results found during the radiation-dominated 
era, this constraint will be satisfied in the infinite past (N = 0) and infinite future (N —i oo) 
for the bouncing solution. Therefore, the matter-dominated solution is expected to describe 
the actual solution accurately only in these two limits. This last point can be easily seen for 
instance for late times. In this case the constraint will be satisfied by a g = 0, which means that 
a = /3af, X = Xf = 1, and thus H = Hf, where 

W 2 = [P^P™ + m4 Al] • (B.2) 

Here we have used eq. (2.6) and taken the N —> oo limit. At late times, when p m /m 4 <C 1, 

the leading order term in this equation will coincide with that of eq. (4.36) and eq. (4.37), but 
the next-to-leading order term will be different as effects of radiation become relevant. During 
early times, the constraint will be satisfied if N = 0, and the same behaviour is found. Only 
the leading order term in the matter-dominated equation accurately describes the evolution of 
the universe. 
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